Open In Colab

01 - Visualizing Neural Responses

Imports

# @markdown Imports
import numpy as np
import matplotlib
matplotlib.rcParams.update({'font.size': 18})
import matplotlib.pyplot as plt
import plotly.graph_objects as go
import numpy as np

# Create figure
fig = go.Figure()

# Add traces, one for each slider step
for step in np.arange(0, 5, 0.1):
    fig.add_trace(
        go.Scatter(
            visible=False,
            line=dict(color="#00CED1", width=6),
            name="𝜈 = " + str(step),
            x=np.arange(0, 10, 0.01),
            y=np.sin(step * np.arange(0, 10, 0.01))))

# Make 10th trace visible
fig.data[10].visible = True

# Create and add slider
steps = []
for i in range(len(fig.data)):
    step = dict(
        method="update",
        args=[{"visible": [False] * len(fig.data)},
              {"title": "Slider switched to step: " + str(i)}],  # layout attribute
    )
    step["args"][0]["visible"][i] = True  # Toggle i'th trace to "visible"
    steps.append(step)

sliders = [dict(
    active=10,
    currentvalue={"prefix": "Bin width:"},
    pad={"t": 50},
    steps=steps
)]

fig.update_layout(
    sliders=sliders
)

fig.show()

Section 1: Experimental set-up

In the next few sections, we will cover important concepts in the context of a case study. In this case study, you are a computational neuroscientist and you have a hypothesis that some neurons in monkey motor cortex correlate with the direction of arm movement while reaching. By this, you mean that a particular neuron will fire more when the monkey reaches in certain directions than others.

You are able to record the spike times of a neuron from motor cortex in a monkey.

This is a typical experimental protocol in neuroscience: recording the spike times of a neuron while presenting the same stimulus over and over, or having the animal perform the same behavior over and over. Each presentation of the stimulus or iteration of the behavior is called a trial. Neurons are noisy - they do not respond exactly the same way every time they see the same stimulus or behavior - collecting multiple trials of data helps us better understand the neural responses and see how variable they are.

You go ahead and perform your experiment and collect a bunch of data. You now have the spike times of the neuron during multiple reachs in the same direction (and you have this for 8 different directions).

Section 2: Raster plots

An excellent and common way to visualize neural spiking in multiple trials is a raster plot, an example of which can be found in Figure 1. In this style of plot, time is represented on the x-axis. Each spike is represented by a vertical line at the time it occurs. Each row (y-axis) can be a separate trial - in this case the raster plot displays a single neuron’s responses over multiple trials and can be used to quickly assess response variability. Note that the timing of each trial needs to be synced. In other words, the stimulus onset, or event onset, should happen at the same time in each trial.

https://github.com/ebatty/EncodingDecodingNotes/blob/main/images/Putnametal.png?raw=True

Figure 1: An example raster plot, adapted from Putnam and Gothard, eNeuro 2019, under CC-BY 4.0

The raster in Figure 1, adapted from Putnam and Gothard, eNeuro 2019, shows neural data from monkey amygdala, recorded while the monkey repeatedly watches a video. The red lines indicate the onset and offset of the video.

Instead of the raster showing one neuron over multiple trials, it can display a population of neuron’s responses, where each row shows the responses of a single neuron. See Figure 3 for an example of this. This figure, from Ito et al, PLoS ONE 2014, represents the spiking responses of simultaneously recorded neurons in cortico-hippocampal brain slices from mice. In this data, there are no trials, this is just recordings over a period of time.

rasterplotexample1

Figure 2: A raster plot with rows as different neurons, from Ito et al, PLoS ONE 2014, under CC-BY 4.0

Section 3: PSTHs

We can summarize the information over trials contained in a raster plot in a peri-stimulus time histogram (PSTH), also sometimes called a peri-event time histogram (PETH). In this plot, we want to show the average underlying firing rate over time. To do that, we need to use discrete chunks of time - we divide time into time bins of some length.

For each time bin, we can then count the number of spikes in that time bin over all the trials. We can then get the average spikes/bin for that time bin during the trial as the the number of spikes divided by the number of trials. We can convert to spikes/second by dividing by the bin length in seconds.

# Define bin width
bin_width = 20

# Create bin edges
bin_edge_times = np.arange(0, trial_length+.1, bin_width)

fig, ax = plt.subplots(1, 1, figsize=(12, 3), sharex=True)

# Make raster of this data
ax.eventplot(all_sp_time_list, colors = 'black')
ax.set(ylabel = 'Trial Number',
      xlim = [0, trial_length],
      ylim = [-.7, 4.7],
      yticks = [0, 1, 2, 3, 4],
      yticklabels = [1, 2, 3, 4, 5]);

# Plot bin vertical lines 
for bin_time in bin_edge_times:
  ax.plot([bin_time, bin_time], [-.7, 4.7], 'r')
../_images/01_VisualizingNeuralResponses_20_0.png

Figure 3: Example raster plot with red vertical lines indicating 20 ms bins

# @markdown 

def plot_raster_and_PSTH(all_sp_time_list, bin_edge_times, PSTH):
  fig, axes = plt.subplots(2, 1, figsize=(12, 7), sharex=True)

  # Make raster of this data
  axes[0].eventplot(all_sp_time_list, colors = 'black')
  axes[0].set(ylabel = 'Trial Number',
        xlim = [0, trial_length],
        ylim = [-.7, 4.7],
        yticks = [0, 1, 2, 3, 4],
        yticklabels = [1, 2, 3, 4, 5]);

  # Plot bin vertical lines 
  for bin_time in bin_edge_times:
    axes[0].plot([bin_time, bin_time], [-.7, 4.7], 'r')

  # Plot PSTH
  x_vals = np.arange(0, trial_length, .1)
  y_vals = np.zeros((len(x_vals)))
  for i_bin in range(len(bin_edge_times) - 1):
    y_vals[(x_vals > bin_edge_times[i_bin]) & (x_vals <= bin_edge_times[i_bin + 1])] = PSTH[i_bin] 

  axes[1].plot(x_vals, y_vals,  'k');
  axes[1].set(xlabel = 'Time (ms)', 
              ylabel = 'Firing rate (spikes/second)');


# Define bin width
bin_width = 20

# Create bin edges
bin_edge_times = np.arange(0, trial_length+.1, bin_width)

# Concatenate all spike times into one array
concat_spikes = np.concatenate(all_sp_time_list, axis=0),

# Get number of spikes in each bin
PSTH, _ = np.histogram(concat_spikes, bins=bin_edge_times) 

# Divide by number of trials and bin length to get spikes/second
PSTH = PSTH/(n_trials * (bin_width/1000))

# Visualize
plot_raster_and_PSTH(all_sp_time_list, bin_edge_times, PSTH)
../_images/01_VisualizingNeuralResponses_26_0.png

Figure 3: PSTH (bottom row) constructed from a raster plot (top row) with a bin width of 5 ms.

# @markdown 

# Compute PSTH
bin_width = 15
bin_edge_times = np.arange(0, 30.1, bin_width)
binned_spikes, _ = np.histogram(np.concatenate(all_sp_time_list, axis=0), bins = bin_edge_times)

# Visualize
plot_raster_and_PSTH(all_sp_time_list, bin_edge_times, binned_spikes)
../_images/01_VisualizingNeuralResponses_30_0.png

Figure 4: PSTH (bottom row) constructed from a raster plot (top row) with a bin width of 15 ms. This bin width is too large and you lose a lot of temporal information!

# @markdown 

# Compute PSTH
bin_width = .5
bin_edge_times = np.arange(0, 30.1, bin_width)
binned_spikes, _ = np.histogram(np.concatenate(all_sp_time_list, axis=0), bins = bin_edge_times)

# Visualize
plot_raster_and_PSTH(all_sp_time_list, bin_edge_times, binned_spikes)
../_images/01_VisualizingNeuralResponses_33_0.png

Figure 5: PSTH (bottom row) constructed from a raster plot (top row) with a bin width of .5 ms. This bin width is too small - the PSTH conveys a lot about the exact spike timings in this set of trials but doesn’t give us a sense of how the underlying firing rate is changing over time

add interactive viz

Let’s return to our proposed experiment. We have recorded the spike times of a motor neuron in a monkey while that monkey reaches in 8 different directions. We have plotted the raster plots given multiple trials of a reach in each direction, see Figure 6 below.

# @markdown

np.random.seed(123)

# Set up parameters for faking data
n_trials = 10
angles = np.array([0, 45, 90, 135, 180, 225, 270, 315])
firing_rates = np.array([5, 8, 40, 54, 55, 38, 24, 8])

# Loop over angles & fake data
reach_sp_times = {}
for i_angle in range(8):
    reach_firing_rates = 15*np.ones((1500))
    reach_firing_rates[(750 - 250):(750 + 250)] = firing_rates[i_angle]

    bin_width = 1/1000

    all_sp_time_list = []
    for i_trial in range(n_trials):

        binned_spikes = np.random.poisson(reach_firing_rates * bin_width)
        sp_times = np.where(binned_spikes)[0]

        all_sp_time_list.append(sp_times)

    reach_sp_times[i_angle] = all_sp_time_list


# Visualize
fig, axes = plt.subplots(2, 4, figsize=(17, 5), sharex=True, sharey = True)
axes = axes.flatten()
for i_angle in range(8):

    axes[i_angle].eventplot(reach_sp_times[i_angle], colors = 'black')
    axes[i_angle].plot([750, 750], [0, 10], 'b')
    axes[i_angle].set(
       ylim = [-.7, 10.7],
       yticks = np.arange(0, 10, 2),
       yticklabels = np.arange(1, 11, 2),
       title = f'{angles[i_angle]} degrees');
    
    if angles[i_angle] > 135:
      axes[i_angle].set(xlabel = 'Time (ms)')
    if angles[i_angle] == 0 or angles[i_angle] == 180:
      axes[i_angle].set(ylabel = 'Trials')

plt.tight_layout()
../_images/01_VisualizingNeuralResponses_38_0.png

Figure 6: Raster plots of motor neuron response for 10 trials for each reach condition. Please note that this is faked data for illustrative purposes. The blue vertical line indicates the start of each reach.